module rk45_mod
    implicit none
    real(kind=8),parameter::kn(5)=[0.25d0, 3.d0/8.d0, 12.d0/13.d0 , 1.d0, 0.5d0]
    real(kind=8),parameter::a2(1)=[0.25d0]
    real(kind=8),parameter::a3(2)=[3.d0/32.d0      , 9.d0/32.d0]
    real(kind=8),parameter::a4(3)=[1932.d0/2197.d0 , -7200.d0/2197.d0 , 7296.d0/2197.d0]
    real(kind=8),parameter::a5(4)=[439.d0/216.d0   , -8.d0            , 3680.d0/513.d0   , -845.d0/4104.d0]
    real(kind=8),parameter::a6(5)=[-8.d0/27.d0     , 2.d0             , -3544.d0/2565.d0 , 1859.d0/4104.d0  , -11.d0/40.d0]
    real(kind=8),parameter::b1(4)=[25.d0/216.d0,1408.d0/2565.d0,2197.d0/4104.d0,-0.2d0]
    real(kind=8),parameter::r1(5)=[16.d0/135.d0, 6656.d0/12825.d0, 28561.d0/56430.d0, -9.d0/50.d0, 2.d0/55.d0]
    real(kind=8) , parameter::eps=1.d-1
    abstract interface
        subroutine func(t,y,yn,n)
            integer,intent(in)::n
            real(kind=8),intent(in)::t
            real(kind=8),intent(in)::y(n)
            real(kind=8),intent(inout)::yn(n)
        end subroutine func
    end interface
contains
    subroutine rk45(f,t,y,n,h,er)
        procedure(func)::f
        integer,intent(in)::n
        real(kind=8),intent(inout)::t
        real(kind=8),intent(inout)::h
        real(kind=8),intent(inout)::y(n)
        real(kind=8)::yn(n)
        real(kind=8)::error(n)
        real(kind=8)::temp(n)
        real(kind=8)::fn(n,6),s,er
        integer::i
        do
            call f(t,y,fn(:,1),n)
            temp=y+a2(1)*h*fn(:,1)
            call f(t+kn(1)*h,temp,fn(:,2),n)
            temp=y
            do i=1,2
                temp=temp+a3(i)*h*fn(:,i)
            end do
            call f(t+kn(2)*h,temp,fn(:,3),n)
            temp=y
            do i=1,3
                temp=temp+a4(i)*h*fn(:,i)
            end do
            call f(t+kn(3)*h,temp,fn(:,4),n)
            temp=y
            do i=1,4
                temp=temp+a5(i)*h*fn(:,i)
            end do
            call f(t+kn(4)*h,temp,fn(:,5),n)
            temp=y
            do i=1,5
                temp=temp+a6(i)*h*fn(:,i)
            end do
            call f(t+kn(5)*h,temp,fn(:,6),n)
            error=y+h*(b1(1)*fn(:,1)+b1(2)*fn(:,3)+b1(3)*fn(:,4)+b1(4)*fn(:,5))
            yn=y+h*(r1(1)*fn(:,1)+r1(2)*fn(:,3)+r1(3)*fn(:,4)+r1(4)*fn(:,5)+r1(5)*fn(:,6))
            er=norm2(yn-error)
            s=((eps*h)/(2*er))**(0.25d0)
            if(s<0.75d0)then
                h=h/2.d0!max(h/2.d0,0.025)
                if(h<1.d-3)then
                    error stop "h < 1.d-3"
                end if
                cycle
            end if
            t=t+h
            y=yn
            if(s>1.5)then
                h=min(2*h,1.d0)
            end if
            exit
        end do
    end subroutine rk45
end module rk45_mod


program main
    use rk45_mod
    implicit none
    real(kind=8)::h
    real(kind=8)::y(3),t,er
    integer::tt
    procedure(func)::f
    y=[0.d0,1.d0,0.d0]
    h=0.1d0
    t=0.d0
    !lorenz curve
    open(10,file="r.txt")
    tt=0
    do
        call rk45(f,t,y,3,h,er)
        write(10,*)y
        write(*,*)t,er
        if(t>10.d0)exit
    end do
    close(10)
end program main

subroutine f(t,y,yn,n)
    integer,intent(in)::n
    real(kind=8),intent(in)::t
    real(kind=8),intent(in)::y(n)
    real(kind=8),intent(inout)::yn(n)
    yn(1)=10.d0*(y(2)-y(1))
    yn(2)=y(1)*(28.d0-y(3))-y(2)
    yn(3)=y(1)*y(2)-8.d0*y(3)/3.d0
end subroutine f
